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^ Abstract 

We analyze the stability properties of the Synchrosqueezing transform, a time-frequency signal analysis method 
that can identify and extract oscillatory components with time-varying frequency and amplitude. We show that 
Synchrosqueezing is robust to bounded perturbations of the signal and to Gaussian white noise. These results justify 
its applicability to noisy or nonuniformly sampled data that is ubiquitous in engineering and the natural sciences. We 

S 



also describe a practical implementation of Synchrosqueezing and provide guidance on tuning its main parameters. 
As a case study in the geosciences, we examine characteristics of a key paleoclimate change in the last 2.5 million 



years, where Synchrosqueezing provides significantly improved insights. 

> 

o 

I. Introduction 

Synchrosqueezing is a time-frequency signal analysis algorithm designed to decompose signals into constituent 
components with time-varying oscillatory characteristics. Such signals f{t) have the general form 

> K 

X f(t) = J2fk(t) + e(t), (1) 

where each component = Ak(t) cos(2tt (f>k(t)) is a Fourier-like oscillatory mode, possibly with time-varying 
amplitude and frequency, and e(t) represents noise or measurement error. The goal is to recover the amplitude 
Ak(t) at the instantaneous frequency (IF) <p' k (t) for each k. 
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Signals of the form ([TJ arise naturally in numerous scientific and engineering applications, where it is often 
important to understand their time-varying spectral properties. Many time-frequency (TF) transforms exist to analyze 
such signals, such as the short-time Fourier transform (STFT), continuous wavelet transform (CWT), and the 
Wigner-Ville distribution (WVD) Q, [11|, [19|, (38j. Synchrosqueezing is related to the class of time-frequency 



reassignment (TFR) algorithms, used in the estimation of IFs from the modulus of a TF representation. TFR 
methods originate from a study of the STFT, which "smears" the energy of the superimposed IFs around their 
center frequencies in the spectrogram. TFR methods apply a post-processing "reassignment" map that focuses the 
spectrogram's energy towards the IF curves and results in a sharpened TF plot. However, standard TFR methods 
do not allow for reconstruction (synthesis) of the components fk(t). (19) , p0| 

Originally introduced in the context of audio signal analysis [15], Synchrosqueezing was recently studied further 



in 1 14 1 and shown to be an alternative to the Empirical Mode Decomposition (EMD) method |24j with a more firm 
theoretical foundation. EMD has been found to be a useful tool for analyzing and decomposing natural signals and, 
like EMD, Synchrosqueezing can extract and delineate components with time-varying spectrum. Furthermore, like 
EMD, and unlike classical TFR techniques, it allows for the reconstruction of these components. Synchrosqueezing 
can be adapted to work "on top of many of the classical TF transforms. In this paper, we focus on the original, 



CWT-based approach studied in [15] and fl4[ , although an STFT-based alternative was developed in [44] and other 
variants are also possible. 

The purpose of this paper is threefold. First, in Section In] we study the stability properties of Synchrosqueezing. 



We build on the theory presented in [14] and prove that Synchrosqueezing is stable under bounded, deterministic 
perturbations in the signal as well as under corruption by Gaussian white noise. This justifies the use of the 
algorithm in real-world cases where different sources of error are present, such as thermal noise incurred from 
signal acquisition or quantization and interpolation errors in processing the data. 



Second, in Section III we explain how Synchrosqueezing is implemented in practice and reformulate the approach 



from 1 14 1 into a discretized form that is more numerically viable and accessible to a wider audience. We also provide 



practical guidelines for choosing several parameters that arise in this process. A MATLAB implementation of the 



algorithm has been developed and is freely available as part of the Synchrosqueezing Toolbox [8J. In Section IV 
we illustrate the algorithm on several numerical test cases. We study its performance and compare it to some of 
the well known TF and TFR techniques. 
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Finally, in Section[Vj we visit a key question in the Earth's climate of the last 2.5 million years (Myr). We analyze 
a calculated solar flux index and paleoclimate records of the oxygen isotope ratio 5 18 0, an index of climate state, 
over this period. We demonstrate that Synchrosqueezing clearly delineates the orbital cycles of the solar radiation 
and provides a greatly improved representation of the projection of orbital signals in <5 ls O records. In comparison 
to previous spectral analyses of 5 18 time series, the Synchrosqueezing representation provides more robust and 
precise estimates in the time-frequency plane, and contributes to our understanding of the link between solar forcing 
and climate response on very long time scales (on the order of lOkyr- lMyr). 

II. The Stability of Synchrosqueezing 

In this section, we state and prove our main theorems on the stability properties of Synchrosqueezing. We first 
review the existing results on wavelet-based Synchrosqueezing and some associated notation and terminology from 



the paper [14]. We define a class of functions (signals) on which the theory is established. 



Definition II.l. [Sums of Intrinsic Mode Type (IMT) Functions] The space A e ,d of superpositions of IMT functions, 
with smoothness e > and separation d > 0, consists of functions having the form fit) = Y^k=ifklf) w i tn 
fk(t) = Akit)e 2m ^ k( ' t \ where for each k, the A k and 4>k satisfy the following conditions. 

A k eL°°nC\ <p k eC 2 , <j/ k , <fi e L°° , inf^(t)>0, 

Vt \A' k (t)\<e\cf>' k (t)\, |$(t)|<e|^(t)|, and 

m + #k-i(t) - ' 

Functions in the class A €t d are composed of several Fourier-like oscillatory components with slowly time-varying 
amplitudes and sufficiently smooth frequencies. The IF components <p' k are strongly separated in the sense that high 
frequency components are spaced exponentially further apart than low frequency ones. 

We normalize the Fourier transform by h(£) = Jf° h(x)e~ 2m ^ x dx and use the notation e = e 1 / 3 . Now for a 
given mother wavelet tp, the continuous wavelet transform (CWT) of / at scale a and time shift b is given by 



W>(o, b) = a' 1 / 2 f{t)^{±=±)dt. If / is supported in (0, oo), then the inversion of the CWT can be expressed 
as f(b) = f °° a~ 3 / 2 Wfia, b)da, where we let = J °° C^iO^ 1 14 P- 6 l We use the CWT to derlne the 
phase transform uif(a,b) by 

dtWf(a,b) n . 
UfM = 2niW f (a,bY (2) 

u)t(a, b) can be thought of as an "FM demodulated" frequency estimate that cancels out the influence of the wavelet 



ij) on Wj(a, b) and results in a modified time-scale representation of /. We can use this to consider the following 
operator. 

Definition II.2. [CWT Synchrosqueezing] Let / 6 A e ^d and h G Cq° be a smooth function such that \\h\\ L1 = 1. 
The Wavelet Synchrosqueezing transform with accuracy 5 and thresholds e and M is defined by 

where Yf~ = {(a,b) : a G [M -1 , M], |W/(a,6)| > e}. We also denote S 5 f -(b,rj) := S 5 f f(b,rj), with the condition 
a G [M _1 ,M] replaced by a > 0. 

For sufficiently small <5, this operator can be thought of as a partial inversion of the CWT of / (over the scale a), 
but only taken over small bands around level curves in the time-scale plane (where Uf(a,b) « 77) and ignoring 
the rest of the plane. As we let 5 — > 0, the domain of the inversion becomes concentrated on the level curves 
{(a, b) : cof(a,b) = r]}. The idea is that this localization process will allow us to recover the components f k more 
accurately than inverting the CWT over the entire time-scale plane. The following theorem was the main result of 



[14]. 



Theorem 11.1. (Daubechies, Lu, Wu) Let f = Ylk=i A k e 27Ti<l}k £ A t) d and e = e 1 / 3 . Pick a function h £ with 
\\h\\ L1 = 1, and pick a wavelet ip 6 C 1 such that its Fourier transform ip is supported in [1 — A, 1 + A] for some 
A < yq^j. Then the following statements hold for each k: 

1) Define the "scale band" Z k = {(a, b) : \a(j)' k {b) — 1| < A}. For each point (a,b) 6 Z k with \Wf(a,b)\ > e, 
we have 

\u f (a,b)-<l/ k (b)\<e, 

and if (a, b) Z k for any k, then \Wf(a, b)\ < e. 

2) There is a constant C\ such that for all b € M, 



lim S s L ,(b, r,)Aj] - A k (b)e 2m ^ 



< Cie. 



-i> J{v-\v-<t>' k (b)\<e} 

This result shows how Synchrosqueezing can identify and extract the components {f k } from /. The first part 
of Theorem 



II. 1 



says that the plot of |>S^-| is concentrated around the instantaneous frequency curves {(f)' k }- The 



second part of Theorem II. 1 tells us that we can reconstruct each component f k by completing the inversion of 



the CWT, locally over small frequency bands surrounding 0', . In particular, it implies that we can recover the 



amplitudes A k by taking absolute values. Theorem II. 1 also suggests that components f k of small magnitude may 



be difficult to detect (as their CWTs become smaller than e). 
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We can now state our new results on the robustness properties of Synchrosqueezing. The following theorem 



shows that the results in Theorem II. 1 essentially still hold if we perturb / by a small (deterministic) error term e. 



Theorem II.2. Let f 6 A 6: d and suppose we have a corresponding e, h, ip and A as given in Theorem II. I 
Suppose that g = f + e, where e £ L°° is a small error term that satisfies \\e\\ Loo < e/ IXLclXf \\ip r i , \\W n ). For 
each k, let M k > 1 be the "maximal frequency range" given by 



M fe =max (^11^11^,(1 + A) 



. Then the following statements hold for each k: 



1 1/2 

1) Assume a 6 [M^ , M k \. For each point (a, b) 6 Z k with \W g (a, b)\ > M k e + e, we have 

\u 9 (a,b) - <p' k (b)\ < C 2 e 

1/2 

for some constant C 2 = 0(M k ). If (a, b) Z k for any k, then \W g [a, b)\ < M fc e + e. 

2) There is a constant C3 = 0(M k ) swc/i that for all 6gM, 

1 



Inn ( ^ / S^ 1/2 (6,f7)d»7 ] - A k (b)e 2m ^ 



< C 3 e, 



Proof: It is clear that 

\W f (a,b)-W g (a,b)\ < ||/-5lLoc a 1 ' 2 



POD 








1 — OO 





(it < a 1/2 e. 



(4) 



Similarly, we also have \di>Wf(a,b) — dt,W g (a,b)\ < a 1//2 e. Now if (a, b) Z k for any k, then using Thm. 
gives 

\W g (a,b)\ < \W g (a,b)-W f (a,b)\ + \W f (a,b)\<Ml /2 e + e. 

1 /2 

On the other hand, if for some k, (a, b) G and | W s (a, 6)| > M fc ' e + e, then by d4i and Thm. 

\u] g (a,b) - <j>' k 



II. 1 



(5) 



II. 1 



< \u g (a,b) - Uf(a,b)\ + \ujf(a,b) - (f>' k (b)\ 

< — . — , — - — ObWf{a, 0) H 



W s (a,&)W/M) 



W fl (a,6) 



+ e 



(^4 

< C 2 e, 



1/2 



e + e e 



H 1- e 



(6) 



where C2 depends only on /, ip and M k . For the second part of Thm. 



n.2 



we fix k and 6 and use the following 
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calculation fl4| p. 12]: 



lim [ Sj ti (b, V )dr,= [ 



(6,/,e,e,oo) 



a - 3/2 W f (a,b)da, 



where 



D(6,/,£i,e 2 ,M) := {a:\W f (a,b)\>e 1 ,\u Jf (a,b)-4>' k (b)\<e 2 ,ae[M- 1 ,M]}. 



(7) 



(8) 



It is also shown in 1 14 p. 12] that if a G D(6, /, e, e, oo), then (a, 6) G Z k , so M fe 1 < a < This means that 



in (jvj), we can replace Sj~(b,r]) by S S j^ Ik (b,ri) and D(6, /, e, e, oo) by D(6, /, e, e, Affc). We can also get a result 
identical to Q for g by simply repeating the argument in 1 14 1. First, note that as 6 — > 0, the expression 



(9) 



'\v-<P' k (b)\<C 2 e 

converges to a~ 3 / 2 W g (a, b)x{\u g (a,b)-<f>' k (b)\<c 2 e}( a ) f° r a l most all a G [M^Mfc], where x is the characteristic 
function of a set. This shows that 



J\V-<t> k (b)\<C 2 e a ' L 



lim 

(a,6)er A/fc *->0J|,,-^(6)|<C 2 6 

g,M ' e + e 



" 3/2 W 5 (a,6)da. 



(10) 
(11) 



'D(fe,<;,M fc 1/2 e+?,C 2 e,M fc ) 

We can justify exchanging the order of integrations and limits in ( [TO] ) by the Fubini and dominated convergence 
theorems, since (9| is bounded by \cT 3 / 2 W g (a, b)\ G L x ({a : \W g (a,b)\ > Adl /2 e + e,a G [M~\M fc ]}) for 
all 5. We also note that |4i and (6 1 show that in the set D(6, /, e, e, M k )\D(b, g, M^ 2 e + e,C 2 e,M k ), we have 
|W/(a,6)| < 2M fc 1/2 f + f. We can now use the result of Thm. II. 1 along with Q, Q and |ll| to find that 



lim 

S->0j\ n -<p> k (b)\<e 



Sf f" (6, rj) drj- lim / 

e ^""^ J\n—d>',. 



\v-^' k (b)\<C 2 e 3,M k e+e 



a - 3/2 W f (a,b) 



< 



D(b,f,e,e,Mk) 

D(6, 9 ,Af t 1/2 £+e-C7 2 e",M t ) 
+ 



D^iW^ e+e,C 2 e,M fc ) 



a _3/2 Wg(a,6)(ia 



- 3 / 2 (W / (a J 6)-W fl (a,6)) 



da 



D(b,f,i,e,M k )\'D(b,g,M 1 k /2 €+e,C 2 e,M k ) 
M k r M k 



- 3 ' 2 W f {a,b) 



da 



< I a~ 1 eda + I a~ J/z [2M^e + ej da 

< (2 log M fc )e + 2 (M fc 1/2 - M~ 1/2 ) (2M fc 1/2 e + e 

< C 3 e. 



-3/2 /W 1 ^ 
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Combining this with the result of Thm. |II.1| finishes the proof. 



Thm. II.2 shows that each component f k can be recovered with an accuracy proportional to the perturbation e and 



its maximal frequency range M k , with mid-range IFs (M& close to 1) resulting in the best estimates. In addition, 



Thm. II.2 implies that we can replace a continuous-time function / with discrete approximations of it. In many 
applications, we only have a collection of samples {f{t n )} available instead of the whole function /, where {t n } 
is a sequence of (possibly nonuniformly spaced) sampling points. We can address this situation in the following 
way. 

Corollary II.3. Let f s G C 2 be the cubic spline interpolant formed from {f(t n )} and define A = sup \t' n+ i — t' n \. 

n 

Then the errors in the estimating <fi' k (b) and fk(b) from f s are both 0(M k A 4: ^ 3 ) for all b. 



Proof: This follows from Thm. II.2 and the following standard estimate on cubic spline approximations [42 
p. 97]: 

ll/.-/l| L -<&A*||/W[U.. 

■ 

This means that we can work with the spline f s instead of /, and as long as the minimum sampling rate A -1 is 
high enough, the results will be close. In practice, we find that the errors are localized in time to areas of low 



sampling rate, low component amplitude, and/or high component frequency (see, e.g., {IV \. 



The second result of this paper is that Sychrosqueezing is also robust to additive Gaussian white noise. We 
start by defining Gaussian white noise in continuous-time. Let S be the Schwartz class of smooth functions with 



rapid decay (see [28]). A (real) stationary generalized Gaussian process G is a random linear functional on S 
such that all finite collections {G(fi)} with /, G S are jointly Gaussian variables and have the same distribution 
for all translates of Such a process is characterized by a mean functional E(G(/i)) = T(/i) and a covariance 



functional E((G(/i) - T(f 1 ))(G(f 2 ) - T(/ 2 ))) = (fi, Rf 2 ) for some operators T : S -> S and R : S -> S, where 



(fi, ia) = /_!° fi(t)f2(t)dt is the 1? inner product. Gaussian white noise N with power a 2 is such a process with 



-oo 

T 2 



T = and R = a I, where / is the identity operator. We refer to |28| for more details on these concepts and to 



[21 1 for basic facts on complex Gaussian variables that are used below. 



Theorem II.4. Let f G A e d and suppose we have a corresponding e, h, if), A and as given in Thm II. I 



and II. 2 with the additional assumptions that tp £ S and < HV'IIl 2 ll^'IU 2 - Let g = f + N, where N is 

Gaussian white noise with spectral density e 2+p for some p > 0. Then the following statements hold for each k: 
1) Assume a £ [Mj~ ,-Mfe]. For each point (a,b) G Z k with \Wf(a, b)\ > e, there are constants E\ and C 2 
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such that with probability 1 — e Eie P , 

\u, g (a,b)-<i/ k (b)\<C 2 e. 

If (a, b) g" Zk for any k, then with probability 1 — e~ E2e F for some constant E 2 , |W fl (a, b)\ < e + \e. 
2) There is a constant C' 3 such that with probability 1 — e~ El<E P , we have for all b Gl that 



2niij> k (b) 



^ W-\v-4>' k {b)\<C' 2 e) 



< C'fe. 



Proof: The CWT of g, W g (a, b), is understood as the Gaussian variable Wf(a, b) + N(ip 0jt i ) ), where ^ a ^{x) = 
a -i/2^ (*=6). We have E(AT(^ a , 6 )) = 0, 

f 2+p poo 



2 +P|U/,l|2 2j 



E(Nty a , b )NU, a , b )) = / V(-)V'(-)^ = e 2+ f(V^)=e 

and since supp(V>) is positive, 

/oo 
-oo 

Similarly, c^W fl (a, 6) is the random variable dfyWf(a, b) + N(tp' a b ), where ^/ b (x) = a~ 3 / 2 ip' (^p). By the same 
arguments as before and noting that supp(-0') C supp(^), we obtain the formulas: 

E(iV« )6 )) = HN«, b ) 2 ) = E(Nty afi )NW afi )) = 

m(< b )m^)) = e 2+p a- 2 \\^\\ 2 L2 

E(iV(^ a , 6 )iV(VQ) = e^a- 1 (^) . 

This shows that the Gaussian variables N(ip a ^) and (N (i/j^b) , N (ip' a b )) G C 2 have zero pseudo-covariance matrices, 
so they are circularly symmetric. If we define the matrix 



V 



then the distribution of (N(ip a ^), N(i()' ab )) is given by 

e -72 i H?("'. 2 )- l/ ~ 1 O. 2 



dwdz. 



7r 2 e 4+2 PdetV 

Since V is invertible and self-adjoint, we can write V^ 1 = U*DU, where D is diagonal and U is unitary. We have 
D 1± D 22 = det(F- 1 ) = det(V)- 1 and D n + D 22 = trace^ 1 ) = (H^H^ + a~ 2 \W\\l 2 ) det(F)- 1 . 
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For a point (a, b), we now define the events G\ = {\N(ip a ^)\ < |}, G 2 = {\N(ij)' ab )\ < f} anc ^ Hk = {l w g( a )^) — 
01(^)1 < f° r eacn ^- We want to estimate P(G\) and P (G\ n G2). Using the above calculations and taking 
E 2 = iHWlt we find that 



2 / re" r dr 



Let Ei = min ag [ A/ -i M ^ g (Z?n + D22) > 0. We note that any rotated polydisk of radius r in (w,z) € C 2 
contains a smaller polydisk of radius 2~ l / 2 r that is aligned with the w and z planes, and use the transformation 

(«/, z') = U(w, z) to estimate 



r e -^{w,z)-V-\w,z) 

p(GinG2) - y w <,, M <„ ^a^v ^ 

/■ e -^P"Kl a +Ai»|*'| a ) 

= / „ 7 2 — — dw'dz' 

J{\(o,i)-u-(w',z')\<^\(i,oyu*(w,z')\<^} 7T e + Pdet 1/ 

/■ e -^(OnK| 2 +fekf) 

> / „ . , „ — - — 77 dw'dz' 

~ J{\ w T+\z'\ 2 <^} vr 2 e 4 + 2 Pdety 

— / — i — 77 dw'dz' 

J{\z'\<2-^e,\w'\<2-^e} vHe 4 +^ det K 

DuD 22 detV J 



re r se s drds 



l_ e -f™ 22 



Now let C 2 = 2M fc 1/2 ||/|| Loo ||V>'|| L i + 3. If (a, 6) Z k for any fc, then Theorem |IU] shows that Gi implies 



|W g (a,6)| < e + |e. Conversely, if (a, 6) G for some fc, we follow the same arguments as in Theorem 
find that 



II.2 



to 



P (H k ) >P {H k \Gi n G 2 ) P (Gi n G 2 



>P 



1 



\W g (a,b)\ 



d b W f (a,b) 



W f (a,b) 



(W g (a, b) - W f (a, b)) - (d b W g (a, b) - d b W f (a, b)) 
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+ e < C' 2 e 



>P 



G 1 r\G 2 jP{G 1 r\G 2 

d b W f (a,b) 



e + ^ + e < G 2 e^ P (Gi n G 2 ) 



>P (2M fc 1/2 ||/|| Loo ||/|| L1 + 3 < G 2 ) P (d n G 2 ) 

=P(G 1 nG 2 ). 



The second statement in Theorem 



II.4 



can be shown in an analogous way. Let C' 3 = 2M^ 2 ((M^ 2 + l)e 2 + 1) + C\ 



and recall the definition ([8]). We fix k and use the above result to estimate 



P 

>P 



5->0 
lim 



lim / |,- <Wk(6 )|<o5gSj;S , i e (6,»?)^-A fc (6)e^W < G^e 

|r?-^)l<^/? " ^ ^ ~ }™ / \v-^(b)\<c^\_ t (b, rj)dq 

H k nG l nG 2 )p(H k nG 1 nG 2 ) 



5->0 
Cie < G^e 



+ 



a- 3/2 N(tl) aib )da + / a^ 3 / 2 VF / (a, 6)da 

D(b, ff ,e-ie,Cie,Af*) JD(b,f,e,e,M k )\D(b,g,e-U,C' 2 e,M k ) 



+ 



Gie < G^e 



tf fe nGinG 2 P(# fc nGinG 2 ) 



Jl4 



3/2^ a + / ' a -3/2 (M l/2 £ + g + + ^ < c , g 



1/Af* 



fffe n Gi n G 2 P(# fc n Gi n G 2 ) 



=P (^2(M fc 1/2 - M" 1/2 )((M fc 1/2 + l)e 2 + 1) + Gi < G3 

=p(GinG 2 ), 



H k n Gi n G 2 P(il fc |Gi n G 2 )P(Gi n G 2 ) 



which completes the proof. 



Part (1) of Theorem II.4 is saying that the noise power gets spread out across the Synchrosqueezing time- 
frequency plane instead of accumulating in a single component instantaneous frequency, despite the fact that the 
Synchrosqueezing frequencies are generally concentrated and are not directly comparable to conventional Fourier 



frequencies (see |44|). Part (2) is the same statement for the entire component A k e lm ^ k , including the amplitude. 
Note that the above argument can also be repeated for more general Gaussian processes such as "1//" noise. In 



this case, the covariances will change (to e.g. ¥,(N(ip a ^)N(t/; a fi)) = e 2+p (ip, Rip)), but the pseudo-covariances 
will still be zero by the translation-invariance of the operator R, and the rest of the argument will be identical. 



11 



III. Implementation Overview 

We now describe the Synchrosqueezing transform in a discretized form that is suitable for efficient numerical 
implementation. We also discuss several issues that arise in this process and how various parameters are to be 



chosen in practice. We are given a vector / G 



2 L+1 , where L is a nonnegative integer. Its elements, 



f m ,m = 0, . . . , n — 1, correspond to a uniform discretization of f(t) taken at the time points t m = to + mAt. To 
prevent boundary effects, we pad / on both sides (using, e.g., reflecting boundary conditions). Figure [j] shows a 
graphical example of each step of the procedure outlined in this section. 




Fig. 1. Synchrosqueezing example for f(t) — (1 + 0.6 cos(2f)) cos(47rf + 1.2t 2 ) and additive noise e(t) ~ A/"(0, 0.5 2 ). Panels in clockwise 
order: a) fit) and f( t) + e (t) sampled, n — 1024 points, b) CWT of /, \Wf\. c) Phase transform uif. d) Synchrosqueezing transform \Tf\; 
with 7 = 10~ 5 (see j ]lII-E| . 



A. DWT of sampled signal: Wj 

We first choose an appropriate mother wavelet if). We pick ift such that its Fourier transform (normalized 



as in Theorem II. li is concentrated in absolute value around some positive frequency £ = uio, and is small and 
rapidly decaying elsewhere (i.e. lim^i^oo P(t)ip(t) = for all polynomials P). Many standard mother wavelets 
satisfy these properties, and we compare several examples in |IV-D 



The DWT samples the CWT Wf at the locations (aj,t m ), where aj = 2^ nv At, j = 1, . . . , Ln v , and the "voice 
number" n v (22j is a user-defined parameter that affects the number of scales we work with (we have found that 
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n v = 32 or 64 works well in practice). The DWT of / can be calculated in 0{n v n\og^n) operations using the 
FFT. We outline the steps below. 

First note that Wf(a, •) = a~ l / 2 i^{— -) ★ /, where * denotes the convolution. In the frequency domain, this 



relationship becomes Wf(a,S t ) = a 1 / 2 f(£,)ip(a£,). We use this to calculate the DWT, Wj(aj,t m ). Let T n (T n 1 ) be 
the standard (inverse) circular Discrete Fourier Transform. Then 



W f {a jr )=T- l [(T n f)®i> j 



(12) 



1/2 ' 

Here denotes elementwise multiplication and ipj is an n-length vector with (ipj) m = a- tp(aj^ m )', Cm are 
samples in the unit frequency interval: £ m = 2irm/n, m = 0, . . . , n — 1. 



B. The phase transform: ujj 

The next step is to calculate the phase transform (|2]). We first require a slight modification of the definition Q, 

1 



"/(a, 6) = ^ r am((W / (a,6))- 1 a 6 W / (a,6)). 



(13) 



In theory Eqs. ( 1_3 1 and ([2]) are equivalent, and in practice (13 1 is a convenient way to obtain a real- valued frequency 



from (|2]). We denote the discretization of ujf by ujj. 

In practice, signals have noise and other artifacts due to, e.g., sampling errors, and computing the phase of Wf is 
unstable when \Wf\ « 0. Therefore, we choose a hard threshold parameter 7 > and disregard any points where 
\Wf\ < 7 . The exact choice of 7 is discussed in Sec. 



III-E 



We use this to define the numerical support of Wj, 



on which ujj can be estimated: 



Sj(m) = [j : \w f ( aj , 



U 



'}• 



> 7 >, for m = 0, . . . , n — 1. 



The derivative in ( [13) can be calculated by taking finite differences of Wf with respect to m, but Fourier 
transforms provide a more accurate alternative. Using the property <%W/(a, £) = 2iri^Wf{a } £), we estimate the 
phase transform, for j G Sj(m), as 

ujj(aj,t m ) = i 3m ^(^/(oj, t m )^) d b Wj(aj,t„ 
with the derivative of W/ estimated via (e.g., (43l) 

W/Cttj-, •) = J"" 1 ((Jn/) © ^) , 

1/2 ^ 

where (dt^j) m = Imaj Crn^(aj^ m )/ At for m = 0, . . . ,n — 1. 

The normalization of 2 corresponds to a dominant, constant frequency of a when /(£) = cos(27rai). This 
allows us to transition from the time-scale plane to a time-frequency plane according to the reassignment map 
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(a, b) — > (w(a, b), b). Note that the phase transform is not the instantaneous frequency itself except in some simple 
cases, but contains requisite "frequency information" that we use to recover the actual frequencies in the next step. 

C. Synchro squeezing in the time -frequency plane: Tf(uj,b) 

We now compute the Synchrosqueezing transform using the reassigned time-frequency plane. Suppose we have 
some "frequency divisions" {wi}^ with wq > and wi + \ > wi for all Let the frequency bin W; be given by 
{w 1 € M : \w' — wi\ < \w' — wi>\ V7' 7^ /}, or in other words, the set of points closer to wi than any other wp. We 
define the discrete-frequency Wavelet Synchrosqueezing transform of / by 

T f (wi,b)= [ W f (a,b)a~ 3/2 da. (14) 



This is essentially the limiting case of the definition ([3]) as 5 — > (note the argument in ([9]) and see also 1 14 p. 



5-6]), but with the frequency variable 77 G R replaced by the discrete intervals Wj. Note that the discretization 
W j is given with respect to n a = Ln v log-scale samples of the scale a, so we correspondingly discretize ( [14] ) 



over a logarithmic scale in a. The transformation a{z) = 2 z / n ", da(z) = a^j^dz, leads to the modified integrand 
W f (a,b)a~ l / 2] ^dz in (JUj). 

To choose the frequency divisions wi, note that the time step At limits the range of frequencies that can be 
estimated. One form of the Nyquist sampling theorem shows that the maximum frequency is w = w na -\ = 
Since / is discretized over an interval of length nAt, the fundamental (minimum) frequency is w = wo = -tj. 
Combining these bounds on a logarithmic scale, we get the divisions wi = 2 lAw w, I = 0, . . . , n a — 1, where 
Aw = -L_l g 2 (n/2). 

We can now calculate a fully discretized estimate of (14i, denoted by TV. Since we have already tabulated uj^ 



and uji[flj,t m ) lands in at most one frequency bin Wi, the integral in ( [14] ) can be computed efficiently by finding 
the associated Wi for each (a,j,t m ) and adding it to the appropriate sum. This results in 0(n a n) computations for 
the entire Synchrosqueezed plane Tj. We summarize this approach in pseudocode in Alg. jlj 

Algorithm 1 Fast calculation of Ty for fixed m 

for I = to n a — 1 do {Initialize T for this to} 

Tf(wi,t m ) <- 
end for 

for all 3 e Sj(m) do {Calculate ([14)} 

{Find frequency bin via wi = 2 lAw w, and cjj e W;} 

I <- min (max (ROUND log 2 ( " /( ^' b "' ) )] , o) , n a - l) 

{Add normalized term to appropriate integral; Az = 1} 

end for 



We remark that as an alternative, the frequency divisions wi can be spaced linearly, instead of the logarithmic 
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scale we use in keeping with the discretization of the CWT in Section III-A| Examples of this approach can be 



found in [45|, but in practice we find that there are no significant differences either way. In principle, the CWT 



itself can be discretized linearly as well, but the approach we took in Section III-A is standard and is preferred for 



its computational efficiency (see |13|, |18|). 



D. Component reconstruction 

We can finally recover each component from Tj by inverting the CWT (integrating) over the frequencies w\ 
that correspond to the kth component, an approach similar to filtering on a conventional TF plot. Let I G £k(t m ) 
be the indices of a small frequency band around the curve of kth component in the phase transform space (based 



on the results of Thm. |II.2| and Thm. |II.4| parts 2). These frequencies can be selected by hand or estimated via a 
standard least-squares ridge extraction method (9J, as done in the Synchrosqueezing Toolbox. Then, using the fact 
that fk is real, we have 

f k (t m ) = ZRj-ffo £ T } {w h t m ) | . (15) 
\leC k (tm) 



where is the normalization constant from Theorem (II. 1 



E. Selecting the threshold 7 

The hard wavelet threshold 7 effectively decides the lowest CWT magnitude at which uj is deemed trustworthy. 
In an ideal setting wherein the signal is not corrupted by noise, this threshold can be set based on the machine 
epsilon (we suggest 1CT 8 for double precision floating point systems). In practice, 7 can be seen as a hard threshold 
on the wavelet representation (shrinking small magnitude coefficients to 0), and its value determines the level of 
filtering. 

In p7| , a nearly minimax optimal procedure was proposed for denoising sufficiently smooth signals corrupted by 
additive white noise. This algorithm consists of soft- or hard-thresholding the wavelet coefficients of the corrupted 
signal, followed by inversion of the filtered wavelet representation. In |T3J, this estimator was also shown to be 
nearly optimal in terms of root mean square error. The asymptotically optimal threshold is y/2 log n ■ a, where 



n is the signal length and a 2 is the noise power. Following |17|, the noise power can be estimated from the 
Median Absolute Deviation (MAD) of the finest level wavelet coefficients. This is the threshold we suggest and 
use throughout our simulations: 

7 = 1.4826^2 logn • MAD(|W^|i :n J 
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where 1.4826 is the multiplicative factor relating the MAD of a Gaussian distribution to its standard deviation, and 
|WV|i :r i„ are the wavelet coefficient magnitudes at the n v finest scales (the first octave). 

IV. Numerical Simulations 

In this section, we provide several numerical examples that illustrate the ideas in Sec. [n] and III and show how 
Synchrosqueezing compares to a variety of other time-frequency transforms in current use. The MATLAB scripts 
used to generate the figures for these examples are available at (8j. 



A. Comparison of Synchrosqueezing with the CWT, STFT and EEMD 

We first compare the Synchrosqueezing time-frequency decomposition to the continuous wavelet transform (CWT) 



and the short-time Fourier transform (STFT) [34]. We show its superior precision, in both time and frequency, at 
identifying the components of complicated oscillatory signals. We then show its ability to reconstruct (via filtering) 
an individual component from a curve in the time-frequency plane. We also compare the recovered component with 



the results of the ensemble empirical mode decomposition (EEMD) method (see [46] for details). 




Fig. 2. Comparison of Synchrosqueezing with the STFT and CWT. (a) Synthetic signal s(t) (bold), corrupted with noise (dashed), shown 
for te [2,8]. (b) STFT of signal s(t). (c) CWT of signal s(t). (d) Synchrosqueezing plot T„(u,t). 
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In Fig. [2] we consider a signal s(t) = si(i) + sz(t) + ss(t) + N(t) denned on t £ [0, 10] that contains different 
kinds of time-varying AM and FM modulation. It is composed of the following components: 

Sl (t) = (l + 0.2cos(t))cos(2vr(2t + 0.3cos(t))), 

s 2 (t) = (l + 0.3cos(2t))e^ /15 cos(2vr(2.4t + 0.5t L2 + 0.3sin(t))) 

s 3 (i) = cos(2vr(5.3t + 0.2t L3 )). 

The signal is discretized to n = 2048 points and corrupted by additive Gaussian white noise N(t) with noise power 
a 2 = 2.4, leading to an SNR of -2.6 dB. 

To make the comparison consistent (as the 7 threshold in Synchrosqueezing has a denoising effect), we first 
denoise the signal using the Wavelet hard-thresholding methodology of § |III-E We then feed this denoised signal to 



the STFT, CWT, and Synchrosqueezing transforms. We use the shifted bump wavelet (see { IV-D[ ) and n v = 32 for 
both the CWT and Synchrosqueezing transforms, and a Hamming window with length 400 and overlap of length 
399 for the STFT. These STFT parameters are selected to have a representation visually balanced between time 
and frequency resolution (34}. 

The component S3 is close to a Fourier harmonic and is clearly identified in the Synchrosqueezing plot T s (Fig. 
Qd)) and the STFT plot (Fig. |2jb)), though the frequency estimate is more precise in T s . The other two components 
have time-varying instantaneous frequencies and can be clearly distinguished in the Synchrosqueezing plot, while 
there is much more smearing and distortion in them in the STFT and CWT. The temporal resolution of the CWT 
and STFT is also significantly lower than for Synchrosqueezing due to the selected parameters. A shorter time 
window or wavelet will provide higher temporal resolution, but lower frequency resolution and more smearing 
between the three components. 

Fig. [3] shows the component S2 reconstructed from the TF plots in Fig. [2] by inverting each transform in a 
small band around the curve of S2- All three time-frequency methods provide comparable results and pick up the 
component reasonably accurately, although the AM behavior around t G [5, 7] is slightly smothered out as a result of 
the noise. On the other hand, EEMD exhibits a poor amplitude recovery and a drifting phase over time, even when 
applied to the original signal without any noise. In general, EMD/EEMD is sensitive to amplitude changes over time 



that impose strong requirements on the frequency separation between the components [37 1, while Synchrosqueezing 



and the other time-frequency methods produce good results as long as the bandwidth of the mother wavelet or 



window is small enough, according to Theorem II. 1 
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Fig. 3. (a-c) Reconstruction of the component S2 on [2, 8] performed by inverting Synchrosqueezing (a), CWT (b) and STFT (c), shown 
as dotted curves, (d-e) The EEMD extraction of S2 with 50 ensembles performed on the signal s (d), and on the same signal without any 
noise and one ensemble (e). The original component S2 is shown in solid curves for reference. 



B. Comparison of Synchrosqueezing with Reassignment Techniques 

We next compare the analysis part of Synchrosqueezing to two of the most common time-frequency reassignment 



(TFR) methods, based on the spectrogram and the Wigner-Ville distribution (see 1 19 1, ch. 4 for details). We apply 
these techniques to s(t), the signal from the last example, for t G [2,8] and with the noise increased to a 2 = 5 
(-5.8 dB SNR). The results are shown in Fig. [4] 

Synchrosqueezing can be understood as a variant of the standard TFR methods. In TFR methods, the directional 
reassignment vector is computed in both time and frequency from the magnitude of the STFT or WVD, which is 
then used to remap the energies in the TF plane of a signal. In contrast, the Synchrosqueezing transform can be 
thought of as a reassignment vector only in the frequency direction. The fact that there are no time shifts in the TF 
plane is what allows the reconstruction of the signal to be possible. We note that in Fig. |4| the Synchrosqueezing 
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Fig. 4. (a) Synchrosqueezing Tj of /. (b) Reassigned spectrogram / STFT of / (RSP). (c) Reassigned smoothed pseudo-WVD of / 
(RWVD). 



TF plot contains fewer spurious components than the other TFR plots. The other TFR methods exhibit additional 
clutter in the TF plane caused by the noise, and the reassigned WVD also contains traces of an extra curve between 



the second and third components, a result of the quadratic cross-terms that are characteristic of the WVD [19|. 



C. Nonuniform Samples and Spline Fitting 

We now demonstrate how Synchrosqueezing analysis and extraction works for a three-component signal that has 
been irregularly sampled. For t G [2,8], let 



f(t) = (1 + 0.5 cos(t)) cos(4vrt) 

+ 2e" 01 ' cos(27r(3t + 0.25 sin(1.4t))) 
+ (1 + 0.5 cos(2.5t)) cos(2vr(5t + 2t L3 )), 



(16) 
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and let the sampling times be perturbations of uniformly spaced times having the form t' m = Atim + At2U m , 
where {u m } is sampled from the uniform distribution on [0, 1]. We take At% = 11/300 and At2 = 11/310, which 
leads to approximately 165 samples on the interval [2, 8] and an average sampling rate of 27.2, or about three 
times the maximum instantaneous frequency of 9.85. As indicated in Cor. H.3| we account for the nonuniform 



sample spacing by fitting a cubic spline through (t' m , f{t' m )) to get the interpolant f s (t), discretized on the finer 
grid t m = mAt with At = 10/1024 and m = 0, . . . , 1023. The resulting vector, f s , is a discretization of the 
original signal plus a spline error term e(t). 

Fig. [HJa-e) shows the Synchrosqueezing TF plot f s and the three reconstructed components. The spline interpolant 
approximates the original signal closely, except for a few oscillations for t > 7.3 where the highest frequencies of / 
occur. The Synchrosqueezing results are largely unaffected by the errors and have no spurious spectral information 
in the TF plot. The effect of the interpolation errors for t > 7.3 is also localized in time and only influences the 
AM recovery of the third component, which contains the highest frequencies and is the most difficult to recover 



as indicated by Thm. II.2 In general, however, we find that components close to the Nyquist frequency are picked 



up fairly accurately as long as the mother wavelet is chosen according to Theorem II. 1 and the components are 



spaced sufficiently far apart (for cases where multiple high frequency components are close together, see the STFT 



Synchrosqueezing approach in [44]). 



D. Invariance to the underlying transform 

As a final example, we show the effect of the underlying mother wavelet on the Synchrosqueezing transform. 



As discussed in |14|, Synchrosqueezing is largely invariant to the choice of the mother wavelet, and the main 
differences one sees in practice are due to the wavelet's relative concentrations in time and frequency (in particular, 
how far away its frequency content is from zero), as opposed to its precise shape. 

Fig. ^ shows the effect of Synchrosqueezing on the discretized spline signal f 8 from the last example, using 
three different complex CWT mother wavelets. These wavelets are: 

a. Morlet (shifted Gaussian) 

MO ocexp(-27r 2 Gu-0 2 ), 

b. Complex Mexican Hat 

MO oc£ 2 exp(-2^V£ 2 ), £>0 

c. Shifted Bump 

MO a exp (-(1 - ((2vr£ - ^/ofy 1 ) , 
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Fig. 5. (a) Nonuniform samples of /, with spline interpolant f s (solid), and original signal / (dashed), (b) Synchrosqueezing TF plane 
TV . (c-e) Extracted components /j? for k = 1,2,3 (solid) compared to originals /& (dashed). 



(e[4-l),4 + l)] 



where for -0 a we use /i = 1, for ijjf, we use a = 1, and for we use fi = 5 and er = 1. These respectively 



correspond to about A = 0.5, 0.25 and 0.16 in Thm. II. 1 We find that, as indicated by Thm. |II. 1] the most accurate 
representation is given by the bump wavelet ip c , whose frequency support is the smallest and exactly (instead of 
approximately) positive and finite. 
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Fig. 6. Wavelet and Synchrosqueezing transforms of / s . Columns (a-c) represent choice of mother wavelet ip a . . . ip c . Top row: |2^>(4f)| 
Center row: |W/ S |. Bottom row: \Tf a \. 



V. Aspects of the mid-Pleistocene transition 

In this section, we apply Synchrosqueezing to analyze the characteristics of a calculated index of the incoming 
solar radiation (insolation) and paleoclimate records of repeated transitions between glacial (cold) and interglacial 
(warm) climates, i.e., ice age cycles, primarily during the Pleistocene epoch (from w 1.8 Myr to « 12kyr before 
the present). The analysis of time series is crucial for paleoclimate research becuase its empirical base consists of 
a growing collection of long deposited records. 



The Earth's climate is a complex, multi-component, nonlinear system with significant stochastic elements [35 |. 
The key external forcing field is the insolation at the top of the atmosphere (TOA). Local insolation has pre- 
dominantly harmonic characteristics in time (diurnal cycle, annual cycle and very long Milankovic orbital cycles) 
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enriched by the solar variability. The response of the planetary climate, which varies at all time scales [26], also 
depends on random perturbations (e.g., volcanism), nonstationary solid boundary conditions (e.g., plate tectonics and 
global ice distribution), internal variability and feedback (e.g., global carbon cycle). Various paleoclimate records, 
called proxies, provide us with information about past climates beyond observational records. These proxies are 
biogeochemical tracers, i.e. molecular or isotopic properties, imprinted into various types of deposits (e.g., deep- 
sea sediment, ice cores, etc.), and they indirectly represent physical conditions (e.g. temperature) at the time of 
deposition. We focus on climate variability during the last 2.5 Myr (that also includes the late Pliocene and the 
Holocene) as recorded by 5 18 The oxygen isotope variations in seawater are expressed as deviations of the ratio of 
ls O to 16 with respect to the present-day standard. Carbonate shells of foraminifera plankton (benthic forams) at 
the bottom of the ocean record 5 18 changes in seawater during their growth. The benthic 5 18 indicates changes 
in the global sea level, ice volume and deep ocean temperature. During the buildup of land ice sheets and the 
decrease in sea level in cold climates, the lighter 16 (9 evaporates more readily than 18 and accumulates in ice 
sheets, leaving the surface water enriched with 18 0. At the same time, the inclusion of ls O during the formation 



of carbonate shells records the ambient seawater temperature of the benthic forams |29|. 



We first examine a calculated element of the daily TOA solar forcing field. Fig. |7Ja) shows f$p, the mid- June 
insolation at 65° N at lkyr intervals J5J. This TOA forcing index has been widely used to gain insight into the 
timing of advances and retreats of ice sheets in the Northern Hemisphere during this period, based on the classic 
Milankovic hypothesis that summer solstice insolation at 65°iV paces ice age cycles (e.g., Q, |23|). The CWT 



and Synchrosqueezing spectral decompositions (using the shifted bump mother wavelet as in the rest of the paper), 
in Fig. [8ja) and Fig. [8je) respectively, show the key time-varying oscillatory components of fgp. Both panels 
confirm the presence of strong precession cycles (at periodicities r=19kyr and 23kyr), obliquity cycles (primary 
at 41 kyr and secondary at 54kyr), and very weak eccentricity cycles (primary periodicities at 95kyr and 124 kyr, 
and secondary at 400 kyr). However, the Synchrosqueezing spectral structure is far more concentrated along the 
frequency (periodicity) direction than the CWT. 

We next analyze the North Altantic and global climate response during the last 2.5 Myr as deposited in benthic 
5 18 in long sediments cores (in which deeper layers contain forams settled further back in time). Fig. [vjb) shows 
fcni- benthic 5 18 0, sampled at irregular time intervals from a single core, DSDP Site 607, in the North Atlantic 
j40j. Fig.gc) shows f C R2- the orbitally tuned benthic 5 ls O stack of (32) (LR05). It is an average of 57 globally 



distributed records placed on a common age model using a graphic correlation technique (3TJ. Fig. |7Jd) shows 



fcR3- the benthic 5 18 stack of [25 1 (H07) calculated from 14 cores (mostly in the Northern Hemisphere) using 
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S(65°N,mid-June) [W/m 2 ] 




-2.5 -2.0 -1.5 -1.0 -0.5 0.0 

Fig. 7. (a) Calculated June 21 TOA insolation flux at 65°7V: fsF- Climate response as recorded by benthic forams S ls O (b) in the 
DSDP607 core /cm, (c) in the LR05 stack f C R2, and (d) in the H07 stack f c R3- 

the extended depth-derived age model free from orbital tuning |27j. The 5 18 records included in these stacks 
vary over different ranges primarily due to different ambient temperatures at different depths of the ocean floor at 
core drill sites. Also, prior to combining the cores in the H07 stack, the record mean between 0.7 Myr ago and the 
present was subtracted from each 6 18 record, so we have different vertical ranges in Fig. |7|b) through Fig. |7jd). 
All 5 18 records are spline interpolated to 1 kyr intervals prior to the spectral analysis. 

The Synchrosqueezing decomposition in the right panels in Fig. [8] is a far more precise time-frequency represen- 
tation of signals from DSDP607 and the stacks than the CWT decomposition in left panels in Fig. [8] or an STFT 
analysis of the H07 stack [25 , Fig. 4]. Noise due to local characteristics and measurement errors of each core is 
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reduced when we shift the spectral analysis from DSDP607 to the stacks, and this is particularly visible in the finer 
scales and higher frequencies. In addition, the stacks in Fig. [HJg) and Fig. [8jh) show far less stochasticity above 
the obliquity band (higher frequencies) compared to DSDP607 in Fig. [8jf). This enables the 23kyr precession 
cycle to appear mostly coherent over the last 1 Myr, especially in comparison to the CWT decompositions. Thanks 
to the stability of Synchrosqueezing, the spectral differences below the obliquity band (lower frequencies) are 
less pronounced betweeen the stacks and DSDP607. Overall, the stacks show less noisy time -periodicity evolution 
than DSDP607 or any other single core due to the averaging of multiple, noisy time series with shared signals. 
The Synchrosqueezing decompositions are much sharper than the corresponding CWT decompositions, and a time 
average of the Synchrosqueezing magnitudes delineates the harmonic components more clearly than a comparable 
Fourier spectrum (not shown). 



During the last 2.5 Myr, the Earth experienced a gradual decrease in the global long-term temperature and 
CO2 concentration, and an increase in mean global ice volume accompanied with glacial-interglacial oscillations 
that have intensified towards the present (shown in Fig. |7Jb) through Fig. |7Jd)). The mid-Pleistocene transition, 
occurring abruptly or gradually sometime between 1.2 Myr and 0.6 Myr ago, was marked by the the shift from 
41 kyr-dominated glacial cycles to lOOkyr-dominated glacial cycles recorded in deep-sea proxies (e.g., Jl2j, |33|, 



[39 1). The cause of the emergence of strong 100 kyr cycle in the late-Pleistocene climate and incoherency of the 
precession band prior to about 1 Myr (evident in Fig. [8jg) and Fig. [8jh)) are still unresolved questions. Both types 
of spectral analyses of selected 5 ls O records indicate that the climate system does not respond linearly to external 
solar forcing. 



The Synchrosqueezing decomposition precisely reveals key modulated signals that rise above the stochastic 
background. The gain (the ratio of the climate response amplitude to insolation forcing amplitude) at a given 
frequency or period, is not constant due to the nonlinearity of the climate system. The 41 kyr obliquity cycle of 
the global climate response is present almost throughout the entire Pleistocene in Fig. [8jg) and Fig. [8pi). The most 
prominent feature of the mid-Plesitocene transition is the initiation of a lower frequency signal («70kyr) about 
1.2 Myr ago that gradually evolves into the dominant 100 kyr component in the late Pleistocene (starting about 
0.6 Myr ago). Finding the exact cause for this transition in the dominant ice age periodicity is beyond the scope of 
our paper, but the Synchrosqueezing analysis of the stacks shows that it is not a direct cause-and-effect response 
to eccentricity variability (very minor variation of the total insolation). 



The precision of the Synchrosqueezing decomposition allows us to achieve a more accurate inversion across 
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any limited frequency band of interest than the CWT spectrum. Inverting the Synchrosqueezing transform over 
the key orbital periodicity bands (i.e. filtering) in Fig. [9] emphasizes the nonlinear relationship between the TOA 
insolation and climate evolution. The top panels in Fig. [9j left to right, show rapidly diminishing contributions 
to the insolation from precession to eccentricity. However, all of the panels below the top row in Fig. [9] show a 
moderately increasing amplitude of variability, i.e., the inverse cascade of climate response with respect to fsF 
from the precession to the eccentricity band in the late-Pleistocene (after 0.6 Myr). On average, the obliquity band 
contains more power than the precession band in DSDP607 and both stacks. Internal feedback mechanisms, most 
likely due to the long-term cooling of the global climate, amplify the response of the eccentricity band after the 
early-Pleistocence (after « 1.2 Myr). The cross-band differences in Fig. [9] (rows 2-4) indicate that a superposition 
of precession cycles can modulate the climate response in lower frequency bands, particularly in the eccentricity 
band, as the climate drifts into a progressively colder and potentially more nonlinear state (e.g., |6[, |36|). 



The Synchrosqueezing analysis of the solar insolation index and benthic 5 18 records makes a new contribution 
in three important ways. First, it produces sharper spectral traces of a complex system's evolution through the 



high-dimensional climate state space than the CWT (e.g., [7, Fig. 2]) or STFT (e.g., 1 12 Fig. 2]). Second, it better 
delineates the effects of noise on specific frequency ranges when comparing a single core to a stack. Low frequency 
components are mostly robust to noise induced by local climate variability, deposition processes and measurement 
techniques. Third, Synchrosqueezing allows for a more accurate reconstruction of the signal components within 
frequency bands of interest than the CWT or STFT. Questions about the key processes governing large-scale 
climate variability over the last 2.5 Myr can be answered by using high-precision data analysis methods such 
as Synchrosqueezing, in combination with a hierarchy of dynamical models at various levels of complexity that 
reproduce the key aspects of the Pliocene-Pleistocene history. The resulting understanding of past climate epochs 



may benefit predictions of the future climate. \41\ 



VI. Conclusions and Future Directions 

Synchrosqueezing can be used to spectrally analyze and decompose a wide variety of signals with high precision 
in time and frequency. An efficient implementation runs in 0(n v n log 2 n) time and is stable against errors in the 
signals, both in theory and in practice. We have shown how it can be used to gain further insight into the climate 
evolution of the past 2.5 million years. 



The authors are also using the Synchrosqueezing transform to study additional topics in climate dynamics, 
meteorology and oceanography (climate variability and change, and large-scale teleconnection), as well as topics 
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in ECG analysis (respiration and T-end detection, fTOj). Synchrosqueezing is also being used by others to address 



problems in the analysis of mechanical transmissions p0| and the design of automated trading systems (TJ. 
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Fig. 8. The CWT time-scale decomposition of (a) the solar forcing index /sf, and the climate response in benthic S 18 of (b) the 
DSDP607 core fcm, (c) the LR05 stack fcn2, and (d) the H07 stack fern- The Synchrosqueezing time-periodicity decomposition of (e) 
the solar forcing index /sf, and the climate response in benthic S la O of (f) the DSDP607 core fcm, (g) the LR05 stack fcn2, and (h) 
the H07 stack fcm- 
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Fig. 9. Milankovic orbital components extracted by inverting the Synchrosqueezing transforms of the insolation index Jsf (a, e, and i), 
and the climate response in benthic (5 18 C> from the DSDP607 core fcm (b, f, and j), the LR05 stack fcR2 (c, g, and k) and the H07 stack 
fcR3 (d, h, and 1). The transforms are inverted over the precession band from 17kyr to 25kyr (left column), the obliquity band from 40kyr 
to 55kyr (middle column), and the eccentricity band from 90kyr to 130 kyr (right column). 



